home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / toolbox / funm.r < prev    next >
Encoding:
Text File  |  1994-09-20  |  1.3 KB  |  53 lines  |  [TEXT/RLAB]
  1. //-------------------------------------------------------------------//
  2.  
  3. //  Syntax:    funm ( A , FUN )
  4.  
  5. //  Description:
  6.  
  7. //  The funm function computes the matrix function of A as described
  8. //  in scalar context by FUN. For example:
  9.  
  10. //    funm ( A , exp )
  11.  
  12. //  calculates the matrix exponential (using expm() is better though).
  13.  
  14. //  The second argument, FUN is NOT a string, it is a variable, in
  15. //  this particular case, it is a variable representing either a
  16. //  builtin, or user-function.
  17.  
  18. //  Funm uses Parlett's method. See Golub and VanLoan (1983), p. 384.
  19.  
  20. //  See Also: expm, logm
  21. //-------------------------------------------------------------------//
  22.  
  23. funm = function ( A , fun )
  24. {
  25.   global (eps)
  26.  
  27.   S = schur (A + zeros (size (A))*1i);    // get Complex Schur form
  28.   F = diag (fun (diag (S.t)));
  29.   tol = eps*norm (S.t,"1");
  30.   n = max (size (A));
  31.   for (p in 1:n-1)
  32.   {
  33.     for (i in 1:n-p)
  34.     {
  35.       j = i+p;
  36.       s = S.t[i;j]*(F[j;j] - F[i;i]);
  37.       if (p > 1)
  38.       {
  39.         k = i+1:j-1;
  40.         s = s + S.t[i;k]*F[k;j] - F[i;k]*S.t[k;j];
  41.       }
  42.       d = S.t[j;j] - S.t[i;i];
  43.       if (abs (d) <= tol)
  44.       {
  45.         fprintf ("stderr", "WARNING: Result from FUNM is probably inaccurate.");
  46.         d = tol;
  47.       }
  48.       F[i;j] = s/d;
  49.     }
  50.   }
  51.   return S.z*F*S.z';
  52. };
  53.